
% this code is designed to work with tracer-filled vessels

function [edgeim, centim] = detectVesselCenters(mov)

numframe = size(mov,3);
edgeim = zeros(size(mov,1),size(mov,2));
centim = edgeim;
counter = 0;
for n = 10:10:numframe
    counter = counter + 1;
    curim = round(mean(mov(:,:,(n-9):n),3));
    edgeim = edgeim + vesselEdges(curim);
    centim = centim + vesselCenterIm(curim);
end
%% normalize both center and edge images:
centim = centim/counter;
edgeim = edgeim/counter;
centim = centim > .4;
edgeim = edgeim > .25;

%% do some morphological operations to clean up the binarized image
centim = bwareaopen(centim,15);
edgeim = bwareaopen(edgeim,15);


centim = bwmorph(centim,'spur',5);
edgeim = bwmorph(edgeim,'spur',5);
%% assign numbers to the vessel segments here:
branchpoints = bwmorph(centim,'branchpoints');
temp = centim~=branchpoints;
temp = bwconncomp(temp);
vesselID = zeros(size(centim));
for n = 1:length(temp.PixelIdxList)
    vesselID(temp.PixelIdxList{n}) = n;
end
centim = bwareaopen(centim,5);
edgeim = bwareaopen(edgeim,5);

centim = centim.*vesselID;      %this stores the vessel IDs of all the center points

function bw = vesselCenterIm(im)

imnorm = im - min(im(:)); imnorm = imnorm/max(imnorm(:));
bw = imbinarize(imnorm, 'adaptive', 'Sensitivity', 0.610000, 'ForegroundPolarity', 'bright');
se = strel('disk',5);
bw = imopen(bw,se);
bw = imfill(bw,'holes');
bw = bwmorph(bw,'thin',inf);
bw = bwmorph(bw,'spur',5);
bw = bwareaopen(bw,10);
return;

function eim = vesselEdges(im)
eim = imgradient(im);
eim = eim-min(eim(:)); eim = eim/max(eim(:));
eim = imbinarize(eim);
eim = bwmorph(eim,'thin',inf);
eim = bwmorph(eim,'spur',4);
eim = bwareaopen(eim,12);
return;